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An efficient pseudo-spectral numerical method is introduced for calculating a self-consistent field (SCF) 
approximation for the linear susceptibility of ordered phases in block copolymer melts (sometimes referred 
to as the random phase approximation). Our method is significantly more efficient than that used in the first 
calculations of this quantity by Shi, Laradji and coworkers, allowing for the study of more strongly segregated 
structures. We have re-examined the stability of several phases of diblock copolymer melts, and find that some 
conclusions of Laradji et al. regarding the stability of the Gyroid phase were the result of insufficient spatial 
resolution. We find that an epitaxial (k = 0) instability of the Gyroid phase with respect to the hexagonal 
phase that was considered previously by Matsen competes extremely closely with an instability that occurs at a 
nonzero crystal wavevector k. 



I. INTRODUCTION 

The local stability of periodic structures, such as those 
formed by block copolymer melts, may be characterized by a 
linear response function that describes the nonlocal response 
of the monomer concentrations to a hypothetical external 
chemical potential field. This response function is closely re- 
lated to the correlation function that is probed by small angle 
x-ray and neutron scattering. The use of a self-consistent field 
(SCF) approximation to calculate the linear susceptibility of 
a homogeneous polymer blend or disordered diblock copoly- 
mer melt is often referred as a "random phase approximation" 
(RPA) for the correlation functionji This usage has become 
entrenched in the polymer literature, despite its somewhat ob- 
scure origin^ '' '* as one of several names for the use of a time- 
dependent SCF theory (i.e., time-dependent Hartree theory) 
to calculate the dynamic linear response of a free electron 
gas^'*. 

In this paper, we present an efficient numerical method to 
calculate the SCF linear susceptibility of ordered phases of 
block copolymers. The SCF susceptibility of the disordered 
phase of a diblock copolymer melt was first calculated by 
Leibler^. Leibler used this both to describe diffuse scatter- 
ing from the disordered phase, and as one building block in 
his theory of weak microphase segregation. Shi, Laradji and 
coworkers**'''''" were the first to calculate the SCF susceptibil- 
ity of ordered phases of diblock copolymer melts, which they 
used to examine the limits of local stability of various ordered 
structures. 

The calculation of the full linear response for an ordered 
structure is a numerically intensive task. Laradji et q/ .^-^'"' 
used a spectral method to complete this calculation that 
was closely analogous to the algorithm used by Matsen and 
Schick" to study equilibrium structures. When applied to 
three dimensionally periodic structures such as the BCC and 
Gyroid phases, this algorithm was able to accurately describe 
only very weakly segregated structures. This limitation was a 
result of a rapid increase in computational cost with increases 
in the number M of spatial degrees of freedom required to 
resolve the structure; The computational cost of a spectral 



calculation of the response to a single perturbation (e.g., the 
response to a single plane wave) is of 0{M^), while the cost 
of a calculation of the full RPA response function (i.e., the 
response to an arbitrary small perturbation) is 0{M'^). 

The earlier use of a spectral method by Matsen and 
Schick" to calculate the equilibrium phase diagram for di- 
block copolymer melts relied heavily upon the use of space 
group symmetry to decrease the number of basis function 
needed to describe a structure. Matsen and Schick introduced 
the use of basis functions with the space group symmetry of 
each structure of interest. Use of these symmetry-adapted ba- 
sis functions functions reduces the number of degrees of free- 
dom M needed to obtain a given spatial resolution by a factor 
roughly equal to the number of point group symmetries in the 
relevant space group. For the BCC and Gyroid cubic phases, 
this reduces M by almost a factor of 48, and thus reduced the 
cost of solution a single iteration of the SCF equations by a 
factor of almost (48)^ « 10^. 

The calculation of the full linear susceptibility, however, 
requires the calculation of the response to arbitrary infinitesi- 
mal perturbations, which generally do not preserve the space 
group symmetry of the crystal. In general, it thus requires the 
use of either a plane wave basis or a spatial discretization that 
does not impose any symmetry upon the perturbation. Pri- 
marily as a result of this loss of the advantages of symmetry, 
Laradji et al.'^'^^ were able to obtain results for the Gyroid 
phase only for xN < 12. Matsen has been able to carry out 
linear response calculations for the BCC and Gyroid phases 
at significantly larger values of xN by considering only the 
response to perturbations that preserve the periodicity of the 
cubic lattice, and that preserve a subgroup of the space group 
of the unperturbed structure'--'^. Matsen's method and results 
are discussed in more detail in Sec. IVIII 

More recent work on numerical methods for solving the 
equilibrium SCFT by Rasmussen and Kalosakas'"^, and by 
Fredrickson and co-worker a'^''^-'^ has made use of pseudo- 
spectral methods for which the cost of a single iteration of the 
SCFT equations scales as 0{M\nM), rather than 0{M^). 
Here, we present a derivation of a perturbation theory for the 
linear response to an arbitrary disturbance in a form that we 
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can evaluate numerically using a pseudo-spectral algorithm. 

The remainder of the paper is organized as follows: In Sec- 
tion [n] we discuss the basic formalism of the SCF linear re- 
sponse theory for perturbations of a periodic microstructure, 
and review the consequences of Bloch's theorem. In Sec- 
tion Hn] we present a perturbation theory for the underlying 
modified diffusion equation, which is used to calculate the 
linear susceptibility of a reference system of non-interacting 
polymers in a periodic potential. In Section HVl we discuss 
the use of linearized SCFT to calculate the corresponding sus- 
ceptibility of an incompressible system of interacting chains. 
In Section [VT] we discuss the implementation and efficiency 
of our algorithm. Section IVIII presents selected results re- 
garding the stability of the HEX, BCC, and (particularly) the 
Gyroid phases of diblock copolymer melts. 

II. RESPONSE OF A PERIODIC STRUCTURE 

In self-consistent field theory, we calculate the average 
concentration pi{r) for monomers of type i by considering 
a hypothetical system of non-interacting polymers in which 
monomers of type i are subjected to a self-consistently de- 
termined field u!i{r). We consider incompressible systems in 
which each monomer occupies a volume v, and define a vol- 
ume fraction field 0j(r) = vpj{r). 

In what follows, we consider the response of an unper- 
turbed state in which oji (r) satisfies the self-consistent field 
(SCF) equation 

w,:(r) =X»j<^j(r)+e(r) . (1) 

Here, Xij is a Flory-Huggins parameter for binary interac- 
tions between monomers of types i and j, and ^(r) is a La- 
grange multiplier field that must be chosen so as to satisfy an 
incompressibility constraint. To calculate the SCF linear sus- 
ceptibility, we consider the perturbation caused by an addi- 
tional infinitesimal external potential Suj°^^{r). The resulting 
deviation 6uji must satisfy 

Scu^ir) = x^JS^J{r) + S^r) + SLu',^\r) (2) 

where (r) is the corresponding deviation in the monomer 
concentration, and (5^(r) is the deviation in ^(r) required to 
satisfy the incompressibility constraint 

Y,^Ur) = (3) 

i 

The deviation S(pi{r) may be described by either of two 
related linear response functions: It may be expressed either 
as an integral 

50,(r) = - J 5„ (r, r')Su;,{r')dr' (4) 

in which Sij (r, r') is nonlocal susceptibility of an inhomoge- 
neous gas of non-interacting polymers to a change in the total 



self-consistent field uj, or by a corresponding relation 

<5</>z = S^J (r, r')<5a;f *(r')dr' (5) 

in which Sij (r, r') is the SCF susceptibility of the interacting 
Uquid of interest to the external potential uj^^*". 

It is convenient to introduce a Fourier representation of the 
problem. As an example, we consider the ideal gas response 
function S in what follows, but identical arguments apply to 
S. The linear response to a perturbation 

-^e^'' '-<5a;,(k) (6) 

k 

may be expressed in Fourier space, for a system of finite vol- 
ume, as a sum 

(50,(k) =^4(k,k')<5c^,(k') (7) 

k' 

in which 

5,, (k,k') = ljdrj dr'5,, (r,r')e-*-+'''' '-' (8) 

where V is the volume of a system containing many unit cells 
of the original structure, with Born-von Karmann boundary 
conditions. 

We are interested here in the response of an unperturbed 
structure that is invariant under translations r r +R, where 
R is any vector in the Bravais lattice of the unperturbed crys- 
tal. As a result, we expect any linear response function to 
exhibit the symmetry 

^„(r,r') = 5,,(r + R,r' + R) (9) 

for any lattice vector R. By Fourier transforming both sides 
of this equality, using definition dHJ for the transform, we find 
that Sij (k, k') can be nonzero only for values of k and k' for 
which e*'^''^'' ) ^ 1 for any lattice vector R. The reciprocal 
lattice is the set of all wavevectors G such that e**^^ = 1 for 
any Bravais lattice vector R. This implies that S'y (k, k') can 
be nonzero only for values of k and k' for which k — k' = G 
for some reciprocal lattice vector G. This is the content of 
Bloch's theorem, as applied to a linear response function. It 
is thus convenient to represent the nonzero elements of S by 
a matrix 

5,j(G,G';k) =5,,(G + k,G' + k) (10) 

where k is a crystal wavevector in the first Brillouin zone. A 
similar notation will be used for the SCF response function 
S. 

Consider the response to perturbation that has the form of 
a Bloch function, 

<5u;,(r) =<5(i,(r)e*'- (11) 
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in which k is a crystal wavevector in the first Brillouin zone, 
and (5a)i(r) is a periodic function with the periodicity of the 
unperturbed lattice, 



in which 



iGr 



G 



(12) 



in which J2g denotes a sum over reciprocal lattice vectors. 
Bloch's theorem guarantees that the resulting density pertur- 
bation will assume the same form 



,(r) =<50,(r)e^' 



(13) 



where S(f){r) is also a periodic function. The relationship be- 
tween the Fourier components ^^(G) and 50i(G) may be 
expressed as a matrix product 



.(G) 



E 

G'j 



5„(G,G';k)(5w,(G') (14) 



with a matrix S'ij(G,G';k) whose elements depend para- 
metrically upon k. It also follows from the block-diagonal 
form of ^(k, k') that the eigenmodes of S and S must be 
Bloch functions. As emphasized by Shi'^, these conclusions 
about the consequences of periodicity are quite general, and 
are independent of the self-consistent field approximation. 

The SCF response function is related to the SCFT free en- 
ergy functional F[(f)] by the standard identity 



S-^iG,G';k) = 



(k + G),50j(-k-G') 



(15) 



Here, the inverse of S is defined in reciprocal space by requir- 
ing 

S;j\G,G';-k)S,k{G' ,G";k) ^ 5,kSG,G" (16) 

The condition for local stability of a periodic structure is thus 
that the matrix S^^^{G, G'; k) be positive definite for every 
k in the first Brillouin zone. The onset of instability oc- 
curs when one of the eigenvalues of ^^^^(G, G'; k) passes 
through zero at some k. 



III. IDEAL GAS RESPONSE 

In SCFT, the monomer concentration fields are obtained by 
calculating the concentrations in a reference system of non- 
interacting chains in which monomers of type i are subjected 
to a field LjJi{r). This reference system is treated by consid- 
ering a pair of constrained partition functions qa{r,s) and 
q|(r, s) for chains of species a, which satisfy the modified 
diffusion equations 



dqg 
ds 

ds 



-Hqa 



(17) 
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(18) 



These quantities satisfy initial conditions q{r, s = 0) = 1 
and ^^(r, s = Na) — 1, respectively, where Na is the length 
of chains of type a. Here, bi and Ui are a statistical segment 
length and a chemical potential field for monomers of the type 
i found at point s along the chains of type a. The volume 
fraction of monomers of type i on chains of type a is given 
by an integral 

cPar{r)^j;^Jdsqa{r,s)ql{r,s) , (19) 

where the integral with respect to s is taken only over those 
blocks that contain monomers of type i. Here, (j)^ is the vol- 
ume fraction of chains of type a, and 



Qa = Jdr qa{r,Na) 



(20) 



The chemical potential fia for species i and Qa are connected 
(for a particular choice of standard state) by a relation 



K = Qae 



(21) 



The SCFT can be applied in either the canonical or grand- 
canonical ensemble by simply regarding either (p^ or /ia as a 
specified input parameter, respectively. 

Here, we consider a perturbation theory for the variation in 
qa that results from a variation in cui. The chemical potential 
field can be written as a sum 



(r , s) = Lul°^ (r , s) + SuJi (r, s) 



(22) 



of an unperturbed part ujf'\r) and a small perturbation 
(5(jJi(r). Similarly <7a(r, s) can be expressed as a sum 



<Za(r, s) =<zi°Hi-, s) + Sqa{r,s) 



(23) 



where gi"'' (r, s) is the solution to the SCF equations for the 
unperturbed periodic structure. Substituting equationsl22]and 
|23]in[T2]yields the perturbation equation: 



— + ) Sqair, s) - -5c^.(r)(?i") (r, s) (24) 



where ij'"' is the unperturbed "Hamiltonian". This must be 
be solved subject to an initial condition Sqa{r,0) = 0. The 
unperturbed fields oj'^^ (r) and gi"-' (r, s) have the periodicity 
and space group symmetry of the unperturbed crystal. 

To take advantage of Bloch's theorem, we consider a 
perturbation of the form of a Bloch function 5uji{r) = 
e^^''^6ijji{r), which will produce a corresponding deviation 



6qair,s) = e'^''6qair,s) 



(25) 
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Here, k is a wavevector in the first Brillouin zone, and SCji 
and 5qa are periodic functions. Substituting these expressions 
into the modified diffusion equation and keeping terms that 
are Hnear in the perturbation yields an inhomogeneous PDE: 



d_ 

ds 



H^]Sqa{r,s) = -SCo^ir)qi°^ir,s) 



in which 



zk)= 



(26) 



(27) 



A closely analogous PDE may be obtained for (5qJ(r, s). 
The resulting deviations must satisfy boundary conditions 
(5qo(r, 0) = and Sql{r, Na) = where Na is the number 
of monomers in a chain of type a. To calculate the ideal gas 
linear response, we numerically solve this pair of PDEs using 
a pseudo-spectral algorithm that is presented in the appendix. 

The perturbation in the periodic part of the monomer con- 
centration field may be expressed, in grand-canonical ensem- 
ble, as an integral 

(5^«(r) = ^ J ^ [%(r,s)gt(r, s) + qair, s)dql{r, s)] 

(28) 

Here, (j)^, Qa, 9a (r, s) and ql^{r, s) all represent values eval- 
uated in the unperturbed state. The integral with respect to s 
in the above must be taken only over the block or blocks that 
contain monomers of type i. 

The expression for (50ai(r) in canonical ensemble is the 
same as that obtained for grand canonical ensemble for any 
crystal wavevector k except k = 0. It may be shown that only 
perturbations with k = (i.e., perturbations with the same 
periodicity as the unperturbed crystal) can induce changes in 
Qa to linear order in the strength of the applied potential. In 
grand-canonical ensemble, any change 5Qa in Qa will cause 
a change 6(t>a ~ SQae^'^/^'^ in the molecular volume fraction 
(pa obtained from Eq. ( |2TI ). but the prefactor to e^" — (pa/Qa 
in Eq. (fT9] l remains constant. In canonical ensemble, where 
4>a is regarded as an fixed input parameter, a change in Qa 
instead induces a change in the denominator of Eq. ( fT9] l for 
(f)ai{v). This yields a slightly modified expression 



.(r) = GCE-0i°)(r)^ 



(29) 



for perturbations in canonical ensemble at exactly k = 0, 
in which "GCE" represents the grand-canonical ensemble re- 
sponse given by the right hand side of Eq. (l28T l. It may 
be shown that this expression yields a perturbation in which 
/ dr 5(i)ai{v) = for all a and i. 

By using the above perturbation theory to calculate the 
Fourier components of the perturbation 54) j (r) caused by a 
particular plane wave perturbation toi(r) oc e**^ we may 
obtain one row of the matrix 5*^^(0, G';k) at a specified 
value of k. The elements of this reciprocal-space matrix are 
generally complex numbers, but may be shown to be real 
when the unperturbed crystal has inversion symmetry. 



IV. SELF-CONSISTENT RESPONSE 

We now discuss how the SCF susceptibility S'(G, G'; k) 
can be calculated from the response function S'(G, G', k) of 
an ideal gas. 

A. General Analysis 

Consider the response to an external perturbation of the 
Bloch form 5Ljf^^{Y) = (5a)™*(r)e''' ''. Substituting the self- 
consistency condition into definition [14] of the ideal gas re- 
sponse function, and expressing the result in Fourier space, 
yields the linear self-consistency condition 



,(G) = -^y(G,G';k)5c:;,(G') 



(30) 



where 



fc,(G') - fcf *(G') + x,k5k{G') + e+S^G') (31) 

Summation over repeated reciprocal wavevectors and 
monomer type indices is implicit. Here we have introduced 
the notation for a vector for which = 1 for all j, and 

5^(G') for a Fourier component of the periodic part S^{r) of 
a deviation 



(32) 



Eq. ( [30l l can be expressed more compactly, in a matrix nota- 
tion, as 



S4>.i = —SijSojj 



(33) 



Here and hereafter, we use boldfaced Greek letters with a sin- 
gle latin monomer index i,j,k, ... to represent column vec- 
tors in the space of reciprocal lattice vectors (or periodic func- 
tions of r), so that d<pi = 6(l)i{G) and Suj = Sijjj{G'), 
and boldfaced capital Roman letter with two monomer type 
indices to represent matrices in this space, so that Sij = 
Sij{G,G';'k). In this notation, matrix-vector and matrix- 
matrix multiplication is thus used to represent summation 
over repeated reciprocal lattice vector arguments. When 
monomer types indices are displayed explicitly, summation 
over repeated indices is implied. 

Imposing the incompressibility constraint 



= J2SMG)=e+SMG) 



yields a condition 



= e+S, 



Solving Eq. ( [35] l for 6$, yields 



(34) 



(35) 



(36) 
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Here, we have introduced the quantities 

5,+ (G,G';k) = 4(G,G';k)e+ 
^+,(G,G';k) = e+5,,(G,G';k) 
5++(G,G';k) = e+5y(G,G';k)e^ 



(37) 



These are represented in matrix notation by Si+, S+j, and 
S++, respectively. Substituting Equation (|36] | for back 
into Equation (l33T l yields 



,50, = -S, 



(38) 



which we use a truncated Fourier representation of M re- 
ciprocal lattice vectors, S is a matrix in a CM dimensional 
space that contains an M dimensional null space (or kernel). 
Thus, though it is tempting for us to rewrite Equation ( 1411 1 as 
= — this would be meaningless, because neither 
S nor S are invertible. The non-null space of the symmetric 
matrix S (i.e., the space spanned by all eigenvectors of S with 
non-zero eigenvalues) is spanned by all functions for which 
= 0' since this is the condition of orthogonality 
with any vector in the null space. The non-null space is thus 
the same as the space of monomer concentration fluctuations 
that respect incompressibility constraint ( [34l i. 



where 



(39) 



The quantity Sij is the SCFT response function of incom- 
pressible system with Xij = 0. 
By solving Eq. dJSl l. we find that 



50, 



where 



i-xs 



kj 



(40) 



(41) 



is the desired SCF response function. Here I denotes the iden- 
tity in the space of reciprocal lattice vectors and monomer 
type indices, with elements SkjScc denotes a matrix in 
this space with elements J^j Xij^jk{G, G'), and inversion is 
defined in this expanded space. 

In order for the incompressibility condition (l34l i to be satis- 
fied for the response to an arbitrary infinitesimal perturbation, 
the response function S for any incompressible liquid must 
satisfy 



= 5+,(G,G';k) = 5,+ (G,G';k) 



(42) 



for any k, G, and G'. The second equality in the above 
follows from Onsager reciprocity. Eq. |42] also implies that 
5++(G, G'; k) = 0. The same conditions apply to S, which 
is a special case of S for x = 0. It is straightforward to con- 
firm that Eq. (|39] l satisfies these conditions for S, and that 
they are preserved by Eq. dTIT i for S. 

Eq. (|42] | implies that S is singular, if viewed as a matrix 
in the space of reciprocal lattice vectors and monomer types, 
since it implies that 



5,,(G,G')V,(G') =0 



(43) 



for any vector of the form ipj{G') = €'^^p{G'), for which 
the elements are functions of G' alone, independent of the 
value of the monomer index j. The matrix Sij thus has a 
null space spanned by the space of all such vectors. In any 
numerical calculation for a system of C monomer types in 



B. Systems with Two Types of Monomer 

Consider a system containing only two types of monomer, 
with a single interaction parameter x — X12 — Xii- In 
this case, it is straightforward to project the problem onto 
the space of physically allowable fluctuations, which respect 
the incompressibility constraint. We define the two compo- 
nent vectors ef = (1, ±1), where i = 1 or 2 for the two 
monomer types, and introduce the following transformation 
for any 2N x 2N dimensional correlation function matrix 
Sij{G, G'): 

5^,(G,G';k) =ef5,,(G,G';k)£^^ , (44) 

in which and i/ belong to the set {+,—}. The quantity 
S'+_i_(G, G') defined in Eq. (l38T l is an example of this no- 
tation. As already noted, the incompressibility constraint re- 



quires that S++ ~ = S I 

remaining element of Eq. 



0. By evaluating the 



for Sij, we find that 



S-+S+i'|_S_|__ 



(45) 



It is straightforward to show that S__ and S__ are related 
by 



= S_ 



I - -S_ 



(46) 



where x = X12 is the conventional scalar Flory-Huggins 
parameter This result is equivalent to Eqs. (17) and (18) 
of Laradji et aZ. '°. The matrix S__ in a system described 
by a truncated basis of M reciprocal lattice vectors and 2 
monomer types is an M x M matrix, which is generally non- 
singular. It therefore is legitimate to rewrite Eq. ( l46b as 



= s: 



(47) 



in close analogy to the RPA equation for an incompressible 
disordered system. 

The limits of stability of an ordered structure may be iden- 
tified by examining the eigenvalues of the RPA response func- 
tion S At each k in the first Brillouin zone, we consider 

the eigenvalue equation 



G' 



(G,G';k)7^„(G';k) = A„(k)V„(G;k) (48) 
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with eigenvectors /„(G; k). The vector ipn{G; k) contains 
the Fourier components of a periodic function V'n(r; k) that 
is the periodic part of an eigenvector of the Bloch form 

e**' '"V'„(r). It follows fromEq. ( |47] i that S and S have 

the same eigenvectors, and that their eigenvalues are related 
by 

A-i(k) = A-i(k)-| (49) 

where A„(k) is a corresponding eigenvalue of the response 
function matrix S The problem of calculating the eigen- 
values of S thus reduces to that of calculating and diagonal- 
izing S. 

C. Response at k = 

When examining hmits of stability we wiU often be partic- 
ularly interested in instabilities at k = 0. In the cases of in- 
terest, these correspond to instabilities toward structures that 
have an epitaxial relationship with the original structure, but 
in which some of the reciprocal lattice vectors of the original 
structure are absent in the final structure. In diblock copoly- 
mer melts, the instabilities of the BCC phase towards hexag- 
onally packed cylinders, and of the cylinder phase towards a 
lamellar phase are found to be epitaxial instabilities of this 
type. The linear response at exactly k = 0, however, has 
some special features that are important to understand when 
constructing a numerical algorithm. 

The response Sij (G, 0; 0) to a perturbation at k = G' = 
corresponds to a response to a spatially homogeneous shift 
in monomer chemical potentials. The change in potential en- 
ergy associated with any spatially homogeneous perturbation 
Su!j^^ in the external fields that couple to monomer concen- 
tration depends only upon the total number of monomers of 
each type in the system, which can change only as a result 
in the change in the number of molecules of each species. 
Such a homogeneous perturbation is thus equivalent to the re- 
sponse to a shift Sfia — J2i ^aiSiLJi in the set of macroscopic 
chemical potential fields for molecules of different species, 
where Nai is the number of i monomers per molecule of 
type a. Such a homogeneous perturbation can have no effect 
upon monomer concentration fields in canonical ensemble, 
because the number of molecules of each type is constrained. 
It also can have no effects in either ensemble in an incom- 
pressible liquid with only one molecular species, such as a 
diblock copolymer melt, because the number of molecules 
per volume is then constrained by incompressibility. 

In either ensemble, we thus expect the matrix 
S--{G, G'; k = 0) in an incompressible diblock copolymer 
melt to have one vanishing eigenvalue, for which the only 
nonzero element of the corresponding eigenvector /„(G; 0) 
is the G = element, corresponding to a homogeneous 
perturbation. This has long been known to be the case in the 
homogeneous phase of an incompressible diblock copolymer 
melt, for which the matrix ^--(G, G'; k = 0) is diagonal. 



with a diagonal element S'(k) = at k = G = G' = 0. 

We have confirmed that our numerical results for both S 

and S__ for ordered phases of diblock melts have this 
property. As in the disordered phase, we also find that the 
eigenvalue A„ (k) associated with one branch of the spectrum 
continuously approaches zero as k ^ 0. To apply Eq. (|49] l 
to a diblock copolymer melt at k = 0, one must thus identify 
and exclude this trivial zero mode. 

The spectrum of S'y (G,G';k — 0) for a three dimen- 
sional structure generally also has three divergent eigenval- 
ues, as a result of translational invariance. The correspond- 
ing eigenvectors are generators of rigid translations, which 
all have the form 

5<^,(r) =5t- V0,(r) (50) 

where 6t is an infinitesimal rigid translation. Basis vec- 
tors for the subspace of rigid translations may be obtained 
considering infinitesimal translations along three orthogonal 
directions. The inverse eigenvalue A^^(k) associated with 
these modes vanish because there is no free energy cost for 
rigid translation of a crystal. In a band structure of A^^(k) 
vs. k for a three dimensional structures, we thus find three 
"phonon-like" bands (one longitudinal and two transverse) in 
which the values of A^^(k) approach zero as fc^ in the hmit 
k^O. 

The fact that these "phonon-like" eigenvalues of Sy di- 
verge as k ^ does not cause any numerical problems for 
diblock copolymer melts if we use Eq. (|49] l to calculate the 

eigenvalues of S In this procedure, we numerically diag- 

onalize the matrix S__, for which the corresponding eigen- 
values have finite values of A,^^(k = 0) = x/2- We find, 
however, that the our numerical results for A~^(k = 0) for 
these modes are very small only when the linear response is 
calculated for a well converged solution to the equilibrium 
SCFT, and only when the calculation is carried out with ad- 
equate spatial resolution. The behavior of these phonon-like 
modes thus provides a useful, and quite stringent, test of nu- 
merical accuracy. 

V. SPACE GROUP SYMMETRY 

When calculating particular eigenvalues of the linear re- 
sponse matrix 5(k) at k = or other special points in the 
Brillouin zone, it is sometimes possible to substantially re- 
duce the cost of the eigenvalue calculation by making use of 
space group symmetry. As a simple example, if we know 
that the instability of a centrosymmetric structure is a result 
of an epitaxial instability towards another centrosymmetric 
structure, we expect the corresponding eigenvector of 5(k) at 
k = to be even under inversion. To calculate the eigenvalue 
associated with such an instability, we may thus calculate a 
matrix representation of S'(k) in the subspace of even func- 
tions by using a basis of cosine functions, rather than plane 
waves, and thereby reduce the number of basis functions by a 
factor of 2 at a given spatial resolution. 
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More generally, group theory can be used in the calcula- 
tion of eigenvalues of iS'(k) at special points in the Brillouin 
zone in a manner very similar to what has long been used in 
the calculation of eigenvalues of the Schroedinger equation in 
band structure calculations. The starting point of the gen- 
eral analysis, in the present context, is the observation that 
the linear response operator S(k) of a periodic structure is 
invariant under the symmetry elements of the so-called "lit- 
tle group" i(k) associated with crystal wavevector k. In the 
case of eigenvectors at k = 0, L(k) is the same as the space 
group G of the unperturbed crystal. More generally, i(k) is 
the subgroup of the full space group G containing symmetry 
elements that leave a plane wave e*'^ '" of wavevector k in- 
variant. By arguments similar to those used to characterize 
the symmetry of eigenvectors of the Schroedinger equation 
at special k points, it may be shown that (in the absence of 
accidental degeneracies) each eigenvalue of S(k) may be as- 
sociated with a specific irreducible representation of group 
L(k). Each such irreducible representation is associated with 
a subspace of functions that transform in a specified way un- 
der the action of the symmetry elements of £(k). Eigenvec- 
tors with the symmetry properties characteristic of a particu- 
lar irreducible representation may be expanded using basis of 
functions that span the associated subspace. 

As a trivial example, consider a one dimensional problem 
involving perturbations at k = of a periodic structure that is 
symmetric under inversion (i.e., a centrosymmetric lamellar 
phase). The space group G of the unperturbed crystal is the 
group —1, which contains only the identity element, r — > r, 
and the inversion element, r — r. At k = 0, the rele- 
vant little group is the same as this full group. There are two 
possible irreducible representations of this group, for which 
the associated subspaces contain all functions ^p{z) that are 
even under inversion, ip{z) = tp{—z), or odd under inversion, 
iplz) ~ —■)/'(— z), respectively. Each eigenvector (z) must 
lie within one of these two subspaces, i.e., must be either even 
or odd. To calculate eigenvalues associated with the subset of 
eigenfunctions at k = that are even under inversion, we can 
use a cosine basis. To obtain eigenvalues associated with the 
remaining odd eigenfunctions, we can use a sine basis. Even 
if we require all of the eigenvalues, the size of the required 
secular matrices is reduced by considering even and odd sub- 
spaces separately, thereby block diagonalizing S. 

In general, to calculate an eigenvector or set of eigen- 
vectors with a known symmetry, we may introduce a 
set of basis function with the desired symmetry. Let 
/i(r), /2(r), . . . , /A/(r) be a set of M orthonormal basis 
functions that lie within the subspace of periodic functions 
associated with a given irreducible representation of the rele- 
vant little group L(k). Each of these basis functions is gen- 
erally a superposition of plane waves with reciprocal lattice 
wavevectors that are related to one another by the symme- 
try elements of group £(k). The phase relationships among 
the coefficients of different plane waves within such a basis 
function are different in different irreducible representations. 
The symmetrized basis functions that we use to represent the 



solution of the unperturbed crystal, which are required to be 
invariant under all elements of the full space group G, are a 
special case of such basis functions, as are cosine and sine 
functions. 

To calculate the response within a subspace spanned by any 
such set of symmetry-adapted basis functions, we consider 
the response of an ideal gas to a perturbation of the form 

<5c^,(r) = e*--^(5^0,/0(r) . (51) 

Such a perturbation is expected to yield a concentration per- 
turbation S(j)a (r) with a periodic factor that can be expanded 
in terms of the same basis functions, with coefficients 6(t)ai- 
The linear response of the ideal gas in the subspace of interest 
may thus be characterized by a matrix 

S(t>ai = -Sal3,ij0i)SuJ/3j (52) 

where Safi,ij (k) is a matrix representation of the ideal gas re- 
sponse S'(k) within the chosen subspace. The corresponding 
RPA response can be represented in this subspace by a matrix 
Sai3,ij0^) of the same form, using the same basis functions. 
The matrix equations that relate the ideal gas and RPA re- 
sponse matrices in this representation are identical to those 
obtained above for the special case of a plane wave basis, ex- 
cept for the replacement of summation over reciprocal vectors 
by summation over basis function indices a and (3. 

Our implementation of the linear response calculation al- 
lows for the introduction of an arbitrary set of such basis 
functions. We have thus far automated the generation of 
symmetry-adapted basis functions, however, only for cases in 
which the required basis functions are invariant under all el- 
ements of i(k), or of a specified subgroup of L(k). In these 
cases, the required basis functions may be generated by the 
same algorithm as that used to generate symmetry adapted 
basis functions for the solution of the unperturbed problem. 
We have thus far actually used symmetry adapted basis func- 
tions, rather than plane waves, only for the purpose of refining 
our results for the eigenvalues of specific eigenvectors of 5* in 
the Gyroid phase, as discussed below. 

When calculating eigenvalues for body-centered crystals 
using a simple cubic cubic computational unit cell, we en- 
countered a subtlety that is a result of the translational sym- 
metry that relates the two equivalent sublattices of the BCC 
structure. This is discussed in appendixlB] 

VI. ALGORITHM AND EFFICIENCY 

To calculate the eigenvalues of S for an equilibrium 
structure of a diblock copolymer melt at a specific crystal 
wavevector k, we first calculate the equilibrium structure, and 
store the converged oj field. If the calculation will use sym- 
metry adapted basis functions, rather than a full plane wave 
basis, we next generate these functions. To calculate the spec- 
trum of S in the invariant space spanned by the chosen set of 
basis functions (which may be plane waves) we must then: 
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1. Use the perturbation theory of Sec. (HIH i to calcu- 
late the ideal-gas response matrix Sij{G, G'; k), in a 
plane wave basis, or Saf3,ij (k), in a basis of symmetry 
adapted functions. 

2. Solve matrix Eq. (l45T l to obtain S 

3. Diagonalize S 

Once the eigenvalues of S__ are known, Eq. (l47l ) may be 

used to obtain the corresponding eigenvalues of S 

In step 1, we calculate Sap^ij using a truncated set of 
M basis function, which may be either plane waves or 
symmetry-adapted functions. To do so, we calculate the per- 
turbation 5(f>i in monomer concentration produced by pertur- 
bations of the form 6ujpj{r) oc e'*' '"/^(r), for every basis 
function fp (r) in our basis set, for perturbations in both uji 
and u!2- This requires us to solve the perturbation theory for 
the ideal gas 2AI times. Each such solution, which requires 
a calculation 6q and 6q^, provides one row of the 2M x 2M 
matrix Sap,ij- 

In the spectral method employed by Shi et al}, 0{M^) 
floating point operations are required to calculate the re- 
sponse of an ideal gas to a single plane wave perturbation. 
The cost of a calculation of the entire matrix 5(G, G'; k) for 
a single k vector is thus O(M^). The matrix operations re- 
quired in steps (2) and (3) each require 0{AI^) operations. 
In the spectral method, the 0{M^) cost of the calculation of 
the ideal gas susceptibility S thus dominates the cost for large 
values of M. With this algorithm, Shi et al— were limited to 
values of M < 800. 

In our pseudo-spectral implementation, the calcula- 
tion of the response to a single perturbation requires 
0{MsMg In Mg) operations, where Mg is the number of grid 
points, or plane waves, and Ms is the number of discretized 
"time-like" steps along the chain. The cost of the calculation 
of the entire matrix ^(G, G'; k) is thus 0{MMsMg InMg). 
For plane wave calculations, M = Mg, and the cost of the 
calculation is 0{MsAP In A/). We have used a time step- 
ping algorithm, described in the appendix, that yields global 
errors of 0{As'^), where As is the contour length step size. 
With this algorithm, very high accuracy can be obtained with 
Ms ~ 10^. The pseudo-spectral algorithm for calculation of 
the ideal gas response function thus becomes much more ef- 
ficient than the spectral method for large values of M. As a 
result, however, the 0{M^) cost of the matrix inversion and 
diagonalization required in steps (2) and (3) will become the 
bottleneck for large M in our algorithm. A comparison of 
CPU times for the different parts of the calculation is given 
in Table |T] The CPU time required for steps 2 and 3 remains 
less than that of step 1 over the range grids reported here, but 
would begin to dominate for slightly larger values of N. 

On a commodity personal computer, the calculation is also 
limited by the memory required in steps 2 and 3. In our imple- 
mentation, in which we have taken care to minimize memory 
usage, these matrix manipulations require storage of 2.5M^ 



double precision real numbers for calculations involving cen- 
trosymmetric unperturbed crystals. 



iV 


M 


time I.G. 


time L.A. 


memory I.G. 


memory L.A. 


8 


512 


2 


0.2 


19 


4 


12 


1728 


30 


7 


35 


55 


16 


4096 


219 


106 


61 


320 


20 


8000 


1053 


771 


107 


1221 



TABLE I: CPU time and memory costs for the calculation of ^(k) 
and its eigenvalues for a Gyroid phase at a given crystal wavevector 
k, using a plane wave basis set on an A^ x x grid, where 
M = A^. Times are in minutes and memory in megabytes. The 
time and memory required to calculate the ideal gas perturbation 
(step 1) are labeled time I.G. and memory I.G., respectively, while 
the costs of the linear algebra operations in steps 2 and 3 are labeled 
time L.A. and memory L.A. Each calculation was carried out on a 
single Athlon 2200 MP processor running at 1.7 GHz. 

The efficiency of our calculation of limits of stability could 
be further improved in several ways that we have not yet ex- 
plored: 

A more efficient algorithm for calculating S (step 2) for 
large values of M might be obtained by replacing our direct 
matrix solution of the linearized SCF equation by an itera- 
tive calculation of the self-consistent field Suj produced by a 
given external field The required iteration would be 

very similar to that which is normally used to solve the SCF 
equations for an equilibrium microstructure. The cost of each 
iteration in such a method would be 0{MsM In M). 

When only a limited number of low eigenvalues are of in- 
terest, as is the case when determining limits of stability, effi- 
cient iterative methods could be used to solve the eigenvalue 
problem. Development of an efficient iterative solution of the 
linear SCF equations would provide a more efficient method 
of calculating the matrix-vector product 6(f> — S6u}°^^ for 
an arbitrary input vector Suj''^^. This could then be used as 
the inner operation of an iterative Krylov subspace method, 
such as the Lanczos method, to efficiently calculate the low- 
est eigenvalues of S. 



VII. STABILITY OF DIBLOCK COPOLYMER PHASES 

In this Section we present "band diagrams" for HEX, BCC 
and Gyroid phases in diblocks. The discussions for HEX and 
BCC phases are focused on the interpretation of the degen- 
erate unstable eigenmodes occurring at k = 0. The dis- 
cussion of Gyroid phase focuses on resolving the questions 
raised by earlier studies of linear stability by Shi, Laradji and 
coworkers^ii^- and by Matsenli. 
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A. Hexagonally Ordered Cylinders 

First, we examine instabilities of the HEX phase toward 
BCC spheres, which occurs at k / 0, and towards a lamellar 
structure, which occurs at k = 0. Experiments have shown 
that diblocks can undergo thermo-reversible transitions be- 
tween cylinders to spheres^'. 

Figure[T]shows the bands of inverse eigenvalues of the SCF 
response function for a HEX phase at / = 0.428 and xN = 
10.9, which lies along the limit of stability of the HEX phase. 
The band structure in this figure appears to agree with that 
obtained by Shi et al. for the same choice of parameters. The 
instability at non-zero k seen in this figure is an instability 
towards a BCC structure, which has been discussed by Shi et 
al. 

Figure |2] shows the evolution of the first few bands in the 
HEX phase with changes in / = 10-9 for a range of com- 
positions near the limit of stability towards a lamellar phase. 
The structure becomes unstable at / = 0.478, when the in- 
verse eigenvalues that are degenerate at k = simultaneously 
pass through zero at k = 0. In addition to these unstable 
modes, this band diagram contains two phonon-like modes, 
for which A^^ (k = 0) = for all choices of parameters. 

Examination of the two unstable eigenvectors at k = 
confirms that they have a structure consistent with an insta- 
bility towards a lamellar phase. Both of the degenerate un- 
stable eigenmodes are found to be even under inversion. We 
find that it is possible to construct three linear superpositions 
of these two eigenmodes such that each has a mirror symme- 
try through one of the three mirror planes of the hexagonal 
phase. 

A more detailed view may be obtained by consider- 
ing the projections of these linear superpositions of the 
unstable modes onto the first "star" of reciprocal vec- 
tors {i.e., the primary scattering peaks) for the HEX 
phase. Let these six primary reciprocal vectors be denoted 
Gi, G2, G3, G4, G5, Gg, with G4 — — Gi, G5 = — G2, 
and Gq = — G3. The projections of the unstable eigen- 
modes onto these vectors can be expressed in real space 
as a sum of three cosine functions with wave-vectors Gi, 
G2, and G3. The amplitudes of these cosine functions for 
the three superpositions discussed above, which each ex- 
hibits a mirror plane, are (1,-0.5,-0.5), (—0.5,1,-0.5), 
and (—0.5, —0.5, 1), respectively. Each of these superposi- 
tions thus tends to increase the amplitude of one of the three 
cosine functions and decrease the amplitudes of the other two 
equally. This is what we expect for an epitaxial instability 
towards a lamellar phase that is aligned along any of three 
equivalent directions. 
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FIG. 1: The bands diagram of inverse eigenvalues of the response 
matrix for the HEX phase. The top subfigure shows the first Bril- 
louin zone and the several points of high symmetry in the hexagonal 
structure^. The lower subfigure shows the band structure of the 
HEX phase at / = 0.428 and = 10.9, with the HEX-BCC 
instability at nonzero k. Subdiagrams labeled T-M and M-K in 
the band diagram show inverse eigenvalues calculated along corre- 
sponding lines in the plane fc^ = in reciprocal space, while subdi- 
agram K-Z shows values along a line constructed perpendicular to 
this plane through point K. 





B. BCC Spheres 

Next, we examine the stability of the BCC spheres and the 
nature of the instability. The instability has been considered 
previously by both Shi et al.^^^ and Matserti?. 



FIG. 2: The lowest few bands for the HEX phase at = 10-9 
and at different compositions near the HEX to lamellar instability, 
for small k along the k-space line segment F-M. This instability 
occurs at k = 0, at / ~ 0.478. 
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FIG. 3: Band diagrams of a stable BCC phase. The top subfigure 
shows the first Brillouin zone and labeled high symmetry points in 
k-space for a BCC crystal.— The lower subfigure shows the BCC 
bands at = 10.7 and / — 0.448, along several line segments in 
reciprocal space that connect the labeled pairs of points T-H, H-P, 
etc. The inverse eigenvalues smaller than 0.15 are plotted. For this 
set of parameters the BCC structure is stable. At a slightly higher 
value of = 10-8 at the same composition the BCC phase be- 
comes unstable, as shown in the Figure|4] 
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FIG. 4: Band diagrams of a BCC phase at xN = 10.8 and / = 
0.448. The inverse eigenvalues smaller than 0.15 are plotted. The 
instability is seen to occur at the F point, i.e., at k = 0. 



The BCC bands for / = 0.448 at two values of xN, 10.7 
and 10.8, calculated using an 8 x 8 x 8 grid, are shown re- 
spectively in Figures [3] and |4] For x^ = 10.7, the BCC 
structure is found stable for all k, whereas as xiV = 10.8, 
the BCC equilibrium structure is unstable. The instability in 
this structure sets in at k = as seen in Figure |4] The unsta- 
ble eigenmode is triply degenerate at k = 0. In addition to 
the unstable mode, the spectrum contains three phonon-like 
bands (one longitudinal and two transverse). The two trans- 
verse phonon bands are degenerate along the lines T-P and 
T-H. 

In this case, we find that it is possible to construct a linear 
superposition of the three degenerate unstable eigenmodes at 
k = such that the resultant superposition has three fold 
symmetry about the [111] axis. The resulting mode has posi- 
tive and equal amplitudes for the 6 primary {011} reciprocal 
lattice vectors that he within a plane perpendicular to per- 
pendicular [111] (thus reinforcing these peaks), and negative 
amplitudes for the remaining 6 {011} vectors (thus leading 
towards their extinction). This superposition coiTesponds to 
a modulation that leads towards the formation of hexagonal 
cylinders along the [111] direction. Equivalent linear super- 
positions can be constructed for instabilities to cylindrical 
phases along the other (111) directions. We thus interpret the 
unstable mode as an epitaxial instability towards a hexagonal 
phase with cylinders along any of the (111) directions. 

This interpretation is consistent with the assumptions un- 
derlying Matsen's calculation, which could only describe in- 
stabilities of this type at k = 0. Our interpretation is, how- 
ever, different from that of Laradji et al. who concluded 
that the instability of the BCC phase was an instability to- 
wards formation of a perforated lamellar structure with lay- 
ers perpendicular to a {011} direction. Laradji et al. did not 
report the fact that this unstable mode is degenerate. It ap- 
pears likely to us that Laradji et al. overlooked the possibility 
of constructing an instability directly to the equilibrium HEX 
phase, rather than a metastable perforated lamellar phase, by 
a suitable linear superposition of the 3 degenerate unstable 
eigenmodes. 



C. Gyroid Phase 

Our results regarding the stability of the Gyroid phase re- 
quire some discussion because earlier results by Laradji et al. 
have been controversial: Laradji, Shi et a/ i^'^''" found that 
Gyroid phase was locally unstable at values of x^ < 12.0, 
within a range of the parameters / and x^ in which the Gy- 
roid phase was then believed to be the equilibrium phase. 
This result, if coiTect, would obviously be incompatible with 
the conclusion that the Gyroid phase was the global minimum 
in the free energy in this range of parameters. Our own inter- 
est in this question was initially raised by the discovery by our 
group22i212^ that an Fddd orthorhombic network is actually 
the equilibrium structure in precisely the slice of parameter 
space (along the line of equal lamellar and HEX free ener- 
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gies) for which Laradji et al. reported the Gyroid to be unsta- 
ble. This raised the question in our minds of whether Laradji 
et al. might have identified an instability of the Gyroid phase 
towards an Fddd phase of lower free energy. (This does not 
appear to be the case, as discussed below). 

Matsen has argued instead that2^ Laradji et al/s conclu- 
sion about the Gyroid phase was probably a result of numeri- 
cal inaccuracy arising from the use of an insufficient number 
of plane waves. Matsen has studied a pathway for epitaxial 
transformations between the Gyroid and a hexagonal cylin- 
der phases in which the cylinders are aligned along the cubic 
[111] axis.— As part of this study, he examined the local sta- 
bility of the Gyroid phase with respect to changes in the com- 
position field that maintained the symmetries shared by the 
Gyroid and cylinder phase: He considered only perturbations 
of the Gyroid phase that retained the periodicity of the parent 
Gyroid phase (corresponding to instabilities at k = 0), in- 
version symmetry, and three fold symmetry around the [111] 
axis. Matsen found that the limit of stability of the Gyroid 
with respect to this type of instability lies near the line of equi- 
librium transitions between HEX and BCC phases, which is 
well beyond the calculated line of equilibrium Gyroid-HEX 
transitions. He thus concluded that the Gyroid phase was lo- 
cally stable with respect to this type of instability throughout, 
and well beyond, the region in which the Gyroid is known to 
have a lower SCF free energy than the HEX phase. 

Matsen assumed throughout this controversy that the in- 
stability found by Laradji et al. must be an instability to- 
wards the HEX phase, of the type that he considered. This 
has remained less clear to us, however, because Laradji et al. 
said nothing about the nature of the instability that they had 
identified, or even whether the instability occurred at zero or 
nonzero crystal wavevector We have confirmed in private 
communications with Shi that he and his coworkers did not 
ascertain the nature of the reported instability. One motiva- 
tion for the work described here was thus to lay this ques- 
tion to rest, by repeating the calculation of the full response 
function without the restrictions on the symmetry or crystal 
wavevector of the perturbation, while using a significantly 
more efficient numerical method than that used by Laradji 
et al.. We find, in agreement with Matsen, that the Gyroid 
phase is locally stable throughout region of parameter space 
in which it has a lower free energy than both the HEX and 
lamellar phases. 

A band diagram for the Gyroid phase at /a — 0.43 and 
XN = 12 is shown in Figure|5] These parameters correspond 
to those at which Laradji et al. concluded that the Gyroid was 
unstable, but that lie within the region in which Matsen and 
Schick found the Gyroid phase to be globally stable (i.e., to 
case 2 of Table 1 in Laradji et aZ. '°). The eigenvalues in this 
diagram were calculated using a 16 x 16 x 16 grid, and a plane 
wave basis. 

The lowest bands in this diagram are phonon-like modes. 
These have a small negative eigenvalue at the T point (k = 0) 
as a result of some remaining numerical error, which is not 
related to the physical instability of the structure. We have 



confirmed that these three eigenmodes at k = correspond 
to rigid translations by confirming that (to within small nu- 
merical errors) they span the same subspace as that spanned 
by Eq. dSOl l for the density modulation generated by arbitrary 
infinitesimal rigid translations. As a corollary of this, they are 
all found to be odd under inversion. We show in Table HIl that 
the associated value of A~^(k — 0) rapidly approaches zero 
as the grid is refined further Because all of the other inverse 
eigenvalues in this diagram, are positive, we conclude that the 
Gyroid is stable at this point in parameter space. 

The physical instability of the Gyroid towards the HEX 
phase that was considered by Matsen is associated with a 
three fold-degenerate eigenvalue at k = that is the next 
band up in this diagram. This set of three degenerate eigen- 
modes is found to have the same symmetry properties (i.e., 
the same irreducible representation) as those found for the 
corresponding instability of the BCC phase towards HEX: All 
three eigenvectors are even under inversion, and the space 
spanned by these eigenvectors contains an eigenvector with 
three-fold rotational symmetry about the [111] axes, as well 
as equivalent eigenvectors with three-fold symmetry about 
each of the other other (111) axes. 

In light of the history of the problem, it is important for 
us to pay attention to questions of numerical convergence. 
Table Hn presents a study of the convergence with increasing 
spatial resolution of the eigenvalues associated with both the 
rigid translation modes (labeled T) and these three dangerous 
modes for the G ^ H instability labeled (H) at k = 0, for 
the same parameters of f — 0.43 and xN — 12. Calculations 
with NxNxN grids of iV=8, 12, 16, and 20 have been car- 
ried out with a plane wave basis. We can only use grids with 
N being a multiple of 4 because the grid must be invariant un- 
der the space group operations of group la'Sd, which include 
diagonal ("d") glide operations that displace the structure by 
one-quarter of a unit cell. Calculations of the eigenvalue of 
the H modes were also carried out on grids with = 16, 
20, 24, and 28 using basis functions with inversion symmetry 
and three fold rotational axis around the [111] axis, which re- 
quire approximately 1/6 as many basis functions at each value 
of N. Calculations of the eigenvalue associated with the rigid 
translation T modes at k = were also carried out for N=16, 
20, 24, and 28 using basis functions with three fold rotational 
symmetry about the [111] axis, with no imposed inversion 
symmetry. This leaves a single rigid translation mode that 
corresponds to a translation parallel to the [111] axis. Calcu- 
lations carried out with a plane wave basis and with symmetry 
adapted basis functions defined on the same grid, for iV = 16 
and N = 20, were found to yield corresponding eigenvalues 
that are identical to within the accuracy displayed in this ta- 
ble. Eigenvalues calculated by different methods are thus not 
distinguished in the table. The values given in the table, and 
those shown in Figs. 1-5, are all values of A~^(k = 0) for a 
diblock copolymer that has been non-dimensionalized by tak- 
ing the reference volume equal to the chain volume, so that 
N = 1. 

For this set of parameters, we conclude that our results are 
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FIG. 5: Band diagrams of a stable Gyroid phase at /=0.43, and 
XN=12, along the representative directions in k-space. The first 16 
inverse eigenvalues are plotted. The labelling of special k- vectors is 
the same as those for BCC structure in Figure |3] 

adequately converged for > 16, with errors of a few times 
lO^'^ for N = 16, but that qualitative errors appear at iV = 8 
and = 12. We find that order of the eigenvalues associated 
with different eigenvectors at k = (which may be uniquely 
identified by their degeneracy and symmetry properties) is in- 
dependent of N for iV > 16, but that the order is different for 
= 8 and N = 12. For iV = 8, at these set of parame- 
ters, we find a total of 12 eigenvectors with negative values of 
A^^(k — 0), in several families of degenerate eigenvectors, 
among which are the 3 degenerate rigid translation/phonon 
modes. At iV = 12, only the 3 translation modes have nega- 
tive eigenvalues. Accurate calculations for the Gyroid phase 
seem to require a substantially finer grid than required for the 
BCC phase, for which we obtain a comparable accuracy at 
similar values of xN using iV = 8. 

The number of plane waves used by Laradji et al. (M < 
783) corresponds most closely to that obtained with an 8 x 8 x 
8 grid (Af — 512), from which we obtain qualitatively incor- 
rect results. The conclusions of Laradji et al. regarding the 
stability of the Gyroid phase at this point in parameter space 
thus appear to be a result of inadequate spatial resolution. 
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12 


16 


20 


24 


28 


T 


-0.2661 


-3.605E-2 


-4.422E-3 


-3.791E-4 


-3.828E-5 


-9.885E-7 


H 


0.1139 


0.1140 


0.1362 


0.1372 


0.1373 


0.1373 



TABLE 11: at k = for the 3-fold degenerate rigid translation 
mode (labelled T) and the 3-fold degenerate dangerous mode for the 
G H instability (labelled H) in the Gyroid phase at / = 0.43 
and xN ~ 12, calculated using different N x N x N grids, for 
iV =8-28. 



We next considered the limits of stability of the Gyroid 
phase with respect to the epitaxial instability that was consid- 
ered previously by Matsen. To do so, we calculated eigen- 



values at k = over a range of values of /a for several 
values of xiV. The three-fold degenerate inverse eigenvalue 
A^^(k = 0) associated with the epitaxial G ^ H insta- 
bility was found to decrease and pass through zero with de- 
creasing J A, at each value of x^, and to be the first inverse 
eigenvalue at k = to become unstable. The resulting limit 
of stability with respect to this eigenmode is shown in Fig- 
ure|6]by the dashed line with open circles. Our results for this 
line of instabilities, which lies very close to the equilibrium 
order-order transition between BCC spheres and hexagonal 
cylinders, agrees very well with Matsen's result for the same 
eigenmodsi^. 

To check whether this instability at k = is preempted by 
another instability at some k 7^ 0, we then calculated another 
band diagram at a point / = 0.3745, x^ = 12 along the pro- 
posed limit of stability line, at which the value of A^^ (k = 0) 
corresponding to the epitaxial G ^ H instability exactly van- 
ishes. The results are shown in Figure|2] At these conditions, 
6 eigenvectors have nearly vanishing eigenvalues at k = 
(the r point), coiTesponding to the 3 rigid translation modes 
and the three G H modes. For this pair of parameters, 
however, we also found very small negative values of A^^ (k) 
at several points near the P and N points: Unstable modes 
were found along the HP(1), PN(2), NT(1) and TP (2) di- 
rections. (The numbers in parentheses indicate the degree of 
degeneracy of the relevant bands along these high-symmetry 
lines.) The most negative eigenmode in this diagram, which 
we assume to correspond to the true limit of stability, lies 
along the HP direction. By calculating eigenvalues at sev- 
eral nearby points along two lines constructed through this 
point along directions perpendicular to the HP line, we have 
confirmed that this is a local minimum of A~^ (k) for the low- 
est band in the three-dimensional band diagram. 

To investigate further, we thus ran a series of calculations 
of the band diagram along the HP line for several of values 
of / near the previously calculated limit of stability of the epi- 
taxial mode, for integer values of xN = 11-16. At each value 
of xN, we found a minimum value of A^^ (k) vs k at approx- 
imately the same wavevector k, which is displaced from the 
P point by a fraction 0.25 of the distance between P and H 
for x^ = 12 and which becomes closer to the H point as 
the value of x^ increases. At each value of xiV, we found 
that this minimum value of A^^(k) at k 7^ passed through 
zero at a value of / that is within approximately 10^'^ of that 
at which the epitaxial k = modes become unstable. For 
example, at x^ = 12, we find an instability at k 7^ at 
/ = 0.3757, vs. / = 0.3745 for the instabihty at k = 0. Dif- 
ferences in the two critical values of / are equally small for 
the other values of xN. We are not confident of our ability to 
accurately resolve such small differences in the critical val- 
ues of / with the 16 x 16 x 16 grid used in these calculations. 
We thus conclude only that this unexpected incommensurate 
instability competes extremely closely with the epitaxial in- 
stability considered by Matsen, and may well preempt it. 

We do not have a physical interpretation of the nature of 
this incommensurate instabihty, or toward what type of struc- 
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FIG. 6: Stability limits of the Gyroid. The full lines are the equi- 
librium phase boundaries (excluding the Fddd phase) and the open 
circles connected by dashed lines are various calculated limits of 
stability of the Gyroid structure. The position of the boundary at 
small values of /, which is an epitaxial (k = 0) instability towards a 
HEX phase, is found to agree well with Matsen's results'^. Another 
eigenmode at k 7^ is found to become unstable along a line that 
would be indistinguishable from this one at the scale of this figure. 
The other symbols o, * and the connecting dashed lines indicate 
the other types of unstable modes, as discussed in the text. 
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FIG. 7: Band diagrams for the Gyroid phase at / = 0.3745, and 
XN — 12, along the representative directions in k-space. The first 
16 inverse eigenvalues are plotted. The labeling of special k- vectors 
is the same as in Figure |5] 

ture it might lead. We note, however, that the identification 
of such an instability at k 7^ could not have been accom- 
plished by Matsen's approach, which allowed him to consider 
only stability with respect to k = eigenmodes with a spec- 
ified residual symmetry: It required the development of an 
efficient method of calculating the full linear response at ar- 
bitrary k. 

In Figure|6]we have also shown three other limits of stabil- 
ity of this Gyroid structure at larger values of /. The branch 



labeled with o is obtained by tracing the equilibrium unit cell 
size a of the Gyroid structure, and is characterized by a rapid 
contraction of the unit cell with increasing /: As this line is 
approached from the left, the derivative da/df appears to di- 
verge to —00. The branch labeled with ★ is obtained from 
perturbation calculation at k = and is doubly degenerate, 
with even eigenfunctions. It is often found to be accompa- 
nied by one of the other two unstable modes not shown in the 
figure, also occurring at k = 0. One of these two modes is 
non-degenerate, for which the eigenfunction is odd under in- 
version, and has three fold rotational symmetry along (111) 
directions. Another one is triply degenerate, with eigenfunc- 
tions that are even. The branch labeled with * is also obtained 
from perturbation calculation at k = and is the result of a 
triply degenerate instability that appears to cross the G H 
instability in the vicinity of the critical point. Along this high- 
/ boundary of the Gyroid phase, where one might expect a 
transition to a lamellar phase, we made no attempt to look for 
yet more instabilities at k ^ 0. 



VIII. CONCLUSIONS 

We have developed a pseudo-spectral algorithm for calcu- 
lating the linear susceptibility of the ordered phases in block 
copolymeric melts that is significantly more efficient than that 
employed in earlier work by Shi and coworkers. We have 
used the new algorithm to re-examine the limits of stability 
of several ordered phases in diblock copolymers, and have 
resolved a controversy regarding the local stability of the Gy- 
roid phase. We have also identified an unexpected instabil- 
ity of the Gyroid phase at a nonzero k vector along the zone 
edge, which competes very closely with the epitaxial G ^ H 
transition considered previously by Matsen. 
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APPENDIX A: INTEGRATION ALGORITHM 

To calculate response of an ideal gas to a specified pertur- 
bation, we must numerically solve Equation ( |26] l for (5g(r, s). 
We do so by discretizing the "time-like" variable s into steps 
of equal size As and numerically integrating the partial dif- 
ferential equation. To carry out the integration, we use a 
pseudo-spectral algorithm closely analogous to one that was 
introduced by Rasmussen and Kalosakas'"* as an algorithm 
for solving the unperturbed MDE that is solved to describe 
unperturbed equilibriums state. We have combined a pseudo- 
spectral algorithm with Richardson extrapolation to obtain 
solutions with errors of C'(As^). 
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1. Unperturbed MDE 

In this subsection, we review the Rasmussen-Kalosakas 
(RK) algorithm for the solution of the unperturbed for q{r, s), 
and present an extrapolation method that we use to improve 
the accuracy of this algorithm. The RK algorithm is based 
upon a representation of ^(r, s) at each value of s on a reg- 
ular spatial grid, and the use of a Fast Fourier Transform 
(FFT) to transform between real- and Fourier-space represen- 
tations. Given a solution q{r,Sn) at s = s„, the value at 
Sn+i — Sn + As may be expressed formally as a product 

q{sn+i) = e-^'"qisn) (Al) 

where q{sn) represents the function g(r, s„), and where 
g-As H jg exponential operator. In the RK algorithm, the 
propagator is approximated by 

The product e^^^''q{sn) is then evaluated by: 

1 . Evaluating a product 

g(+)(r) = e-^-Wg(r,s„) (A3) 
at regularly spaced grid points, 

2. Applying a FFT to obtain g+ (k) and using the Fourier 
representation to evaluate 

gi;\(k)=e-^^*l''l^(+)(k) (A4) 

3. Applying an inverse FFT to obtain q^^\(r) and again 
evaluating a product 

(?(r,s„+0 = e-^-Wq(-\(r) , (A5) 
on the real-space grid. 

This algorithm yields a solution with local errors of 0{As^), 
or global errors of O(As^) 

To improve the accuracy of the solution we have used an 
extrapolation scheme in which we calculate each time step 
using two different values of As, which differ by a factor of 
2, and then extrapolating to As = to obtain the next value. 
Given g(s„), we first calculate a function (7(s„+i; As) by ap- 
plying the RK algorithm once with a time step As. We then 
calculate a function q(s„+i; As/2) by applying the above al- 
gorithm twice, using a step size As/2. The final value of 
q{sn+i), which is used as the starting point for the next step, 
is obtained from the extrapolation 

g(s„+i) = [4g(s„+i; As/2) -q(s„+i; As) ]/3 , (A6) 

which is designed to cancel the accumulating errors of order 
(As)2. 



This extrapolation scheme for the unperturbed MDE was 
originally implemented in the pseudo-spectral version of our 
SCFT code. When designed, it was expected to remove 
global errors of order ©(As^), and to leave errors of ©(As^). 
When the algorithm was tested, however, it was found to 
yield a global error that decreased with As as (As)"'. We 
now believe that this behavior is a result of a special prop- 
erty of reversible integration algorithms^'. An algorithm is 
said to be reversible if one forward propagation of the so- 
lution q{sn), followed by a backward propagation with the 
same step size, can exactly recover the starting point q{sn)- 
The RK algorithm is reversible in this sense. It is known^- 
that the Taylor series expansion of the global error produced 
by a reversible discrete integrator for any system of first or- 
der differential equations contains only even powers of As. 
As a result, reversible integration algorithms always exhibit 
the behavior we observed for the extrapolated RK algorithm: 
An extrapolation that is designed to decrease the global errors 
from ©(As)^" to order (As)^"+'^ generally yields a solution 
with global errors of order (As)^"+^. 



2. Perturbed MDE 

We now describe the algorithm used to solve Eq. ( |26] |. 
and the conjugate equation for Sq'^ . A formal solution of the 
Equation|26]over a single step of length As can be written as 
an integral: 

5g(s„+i) - Gk(As)5g(s„) + / ds'Gk(s')/(s„+i - s') 

Jo 

(A7) 

where 

/(r,s) = -fc(r;k)qo(r,s) (A8) 

is the inhomogeneous term in this linear PDF, and where 
G'k(s) = e^'^^^ is the propagator for perturbations of crystal 
wavevector k. 

Our integration scheme (prior to Richardson extrapolation) 
is to take 

6q{s„+i) = G{As)dq{sn) + G(As/2)7(s„)As (A9) 
where 

7(S„) EE [/(s„+i)+/(s„)]/2 (AlO) 

while using the RK approximation for propagation by both 
G(As) and G(As/2). This algorithm, like the RK algorithm, 
yields global errors of 0{As^). Also, like the RK algorithm, 
it is reversible. We thus use the same extrapolation scheme for 
this algorithm as that given in Eq. ( IA6b for the RK algorithm. 
After 6q and 6q^ are obtained by this method, the integral 
with respect to s in Eq. (|28] | is evaluated using Simpson's 
method. The extrapolated solution for S4>, and thus for the 
matrix S, have errors of 0{As^), as demonstrated in Fig. |8] 
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FIG. 8: First inverse eigenvalue of the response function S a lamellar 
phase vs. As* when the perturbation theory is evaluated using the 
extrapolation scheme of Eq. IA6I 



BCC structure using a simple cubic cell, without explicitly 
accounting for the centering symmetry of the unperturbed 
structure, the list of eigenvalues obtained at a specified k 
in the FBZ of the simple cubic cell include both those as- 
sociated with k and those associated with another wavevec- 
tor k' = k + G that differs from k by a lattice vector 
G = (±1, 0, 0)27r/a, (0, ±1, 0)27r/a, or (0, 0, ±l)27r/a that 
is part of the simple cubic reciprocal lattice, but not part of 
the reciprocal lattice of BCC. The vectors k and k' are thus 
equivalent from the point of view of the simple cubic lattice 
but inequivalent from the point of view of the BCC lattice. 
Correspondingly, if k lies in the FBZ of the simple cubic cell, 
then k' generally lies within the FBZ of the BCC unit cell but 
the outside the smaller FBZ of the simple cubic unit cell. The 
result is a "folding" of the FBZ of the BCC unit cell into the 
smaller FBZ of the simple cubic computational unit cell. 



APPENDIX B: BODY- AND FACE-CENTERED LATTICES 

Here, we explain an issue that we encountered when cal- 
culating band diagrams for BCC and Gyroid phases using 
a non-primitive cubic unit cell. Analogous issues can arise 
whenever a body- or face-centered crystal is treated with a 
non-primitive computational unit cell. 

Band diagrams for the BCC and Gyroid phases, can be cal- 
culated using either a cubic unit cell with orthogonal axes, or 
a primitive unit cell with non-orthogonal axes, which has half 
the volume of the cubic unit cell. (The Gyroid structure is 
based on a BCC lattice). We have used a simple cubic com- 
putational unit cell with cell size a, and discretized this with 
a simple cubic FFT grid. The reciprocal lattice associated 
with this computational unit cell thus includes wavevectors 
that are not part of the reciprocal lattice of the BCC Bravais 
lattice, which is an FCC lattice in fc-space. The first Bril- 
louin zone (FBZ) associated with the simple cubic unit cell 
(|fca:|, |fcy|, |fcz| < 7r/2) has half the volume in /c-space as the 
FBZ for the BCC unit cell. 

If the algorithm described in this paper is applied to a 



The problem could be avoided either by using a primi- 
tive BCC unit cell from the outset, or by using only the re- 
ciprocal vectors of the BCC unit cell in the calculation of 
Sij (G, G'; k). What we have actually done is to use our ma- 
chinery for generating symmetrized basis functions to gener- 
ate basis functions that are invariant under a subgroup of L{k) 
group that, at a minimum, includes the identity r ^ r and the 
body-centering translational symmetry r ^ r + (1, 1, l)a/2. 
This guarantees that we will only obtain eigenfunctions of 
the Bloch form e**' '"V'Ti(r) in which V'Ti(r) has the period- 
icity of the BCC lattice, and not just the periodicity of the 
larger cubic cell. If only these two symmetry elements are in- 
cluded, the algorithm automatically generates basis functions 
that are simply plane waves with wavevectors that belong to 
the reciprocal lattice for BCC, while discarding the remaining 
"extinct" reciprocal lattice vectors of the simple cubic lattice. 
The advantage of this approach is its generality: It allows the 
use of either primitive or non-primitive unit cells, generalizes 
immediately to the treatment of other kinds of body- and face- 
centered crystals, and does not require the explicit addition of 
special extinction conditions to our algorithm. 
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